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Four generalizations of the Phase Integral Approximation (PIA) to sets of ordinary differential 
equations of Schrodinger type {uj{x) + '^'^^-^Rjk{x)uk{x) = 0, j — 1,2, ... are described. 
The recurrence relations for higher order corrections are given in a form valid to arbitrary order and 
for the matrix R(a;)(= {Rjkix)}) either hermitian or non-hermitian. For hermitian and negative 
definite R(x) matrices, a Wronskian conserving PIA theory is formulated which generalizes FuUing's 
current conserving theory pertinent to positive definite R(a::) matrices. The idea of a modification 
of the PIA, well known for one equation {u" {x) + R{x)u{x) = O) is generalized to sets. A sim- 
plification of Wronskian or current conserving theories is proposed which in each order eliminates 
one integration from the formulas for higher order corrections. If the PIA is generated by a non- 
degenerate eigenvalue of the R(2;) matrix, the eliminated integration is the only one present. In 
that case, the simplified theory becomes fully algorithmic and is generalized to non-hermitian R(2;) 
matrices. General theory is illustrated by a few examples generated automatically by using the 
author's program in Mathematica published in larXiv:0710. 54061 [math-ph]. 



I. INTRODUCTION 



This paper deals with generalizations of the well known Phase Integral Approximation i^iii^i^il i I2ii3,i6,i7 fpjj^is ^p. 
proximation was developed for solutions of the one-dimensional time independent wave equation, 

u"{x)+R{x)u{x)=0, (1) 

(e.g., R{x) — j^[E ~ V{x)] for the Schrodinger equation in cartesian coordinates). Possible generalizations of this 
theory to sets of ODEs of similar type: 

N 

u';{x)+Y.^Jf'-^^>k{x)^0, j^l,2,...,N, (2) 
fc=i 

will be described. This can be regarded as going from a "scalar case", Eq. ([1]), to a "vector case": 

u"(x) -I-R(x) • u(x) = 0, (3) 

where vector \i{x) and matrix R(x) have elements Uj{x) and Rjk{x), j, fc = 1, 2, . . . , and the dot (here and in what 
follows) denotes summation over neighboring indices of vectors and/or matrices (contraction). Note also that in this 
paper the convention is adopted that the prime indicates differentiation with respect to the variable indicated in the 
argument of the function in question. 

Basic results of scalar theory in a form convenient for generalizations are described in Sec. [TTl This theory was first 
generalized to vector cases with a hermitian positive definite R matrix in 1979 by S. A. Fulling.* In Sees. IIIII and llVi 
FuUing's results will be presented in a somewhat modified and more general form and extended to negative definite 
R(a;) matrices. The original treatment will also be commented on briefly. Furthermore, a simplified PIA theory will 
be proposed, containing no integrals characteristic for the current or Wronskian conserving theories. 

In lowest order, the phase integral approximation for a two dimensional vector case {N = 2) was also introduced 
independentlyi^ and then generalizedi^ to any A^ > 1. This theory was useful in providing initial conditions for 
numerical integration of the two relevant differential equations. The eigenvalue problem for these ODEs was solved 
numerically in the limit in which the numerical integration interval tended to infinity and the accuracy required was 
very high. This calculation would not be possible without efficient asymptotics at x — s- cx) provided by the phase 
integral approximation. Extension of this theory (valid also if R(a;) is non hermitian) to higher orders is possible but 
turns out to be rather complicated. A simpler theory of possibly non hermitian vector cases is given in Sec. |Vl 

In Sec IVI[ for the simplest vector case of A^ = 2, all earlier discussed vector theories are compared. 

Sec. IVIII describes singularities in the PIA and a possible choice of the auxiliary function a{x) present in all earlier 
theories. 

Sec. IVIIll gives examples produced by the author's program in Mathematical" and Sec. IIXI contains conclusions. 
In the rest of this introductory section, we discuss a few simple facts pertaining to exact solutions of Eqs. (H)) and 
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Equations of the form ^ can be arrived at from somewhat more general "Schrodinger hke" equations: 



N 



1'^ (x) + aj {x)u'j (x) + J2 Rjk {x) Uk{x)=0. (4) 



k=l 



Uj(a;) = exp ^ / aj{x)dx Uj{x), (5) 



Using the transformation: 



the first derivative terms are ehminated, and equations ([¥]) take the form ^ with 

R,k{x) = R.kix) - 5,k \ [\a%x) + a'^ix)]. (6) 

For the radial part of the Schrodinger equations in spherical coordinates, x — r (the spherical radius), and aj{r) — 2/r, 
leading to 

u,{r) = ru,{r), R,,{r) ^ R,,{r), R^.^[r) ^'^-^[E -V{r)] (7) 

If the function R{x) or the matrix R(x) is real for real x, as often happens in applications, we can assume that the 
solution u{x) or u{x) is also real. Dealing with complex solutions is either a question of convenience (e.g., complex 
exponential solutions for constant R{x) or R(a;)) or is due to the physical meaning of the solution (e.g., a wave 
function in quantum mechanics). In any case, however, the real or imaginary part of a complex solution u{x) or u(x) 
is also a solution. 

In the vector case, an important special situation arises if the matrix 'R,(x) is hermitian ("hermitian vector case"): 

R,k{x) ^ Rl^{x). (8) 
In that case the eigenvalues of R(a::) are real, and the exact solution of Eq. ^ conserves the generalized current>S 

N 

iTuT — n. rrnT = Tm 

dx 



(Tat = 0, a]\j = lm'^^u*{x)u'j{x) =lm[u{x),u'{x)), (9) 

where the compact notation is obtained if one introduces the scalar product in the N dimensional complex Hilbert 
space 7i^, (a, b) = a* - b = X^^Li '^j^j- For = 1, cti is proportional to the quantum mechanical probability current 

An important subgroup of hermitian vector cases is the situation where R(x) is both hermitian and real ("real 
hermitian case" ) , 

linRjk{x) = and Rjk{x) = Rkj{x). (10) 

In that case, the eigenvalues of R(x) are again real but furthermore we can also assume that the eigenvectors of R(a;) 
are real. 

Note that both the hermitian and real hermitian vector cases can be considered to be generealizations of those 
important scalar cases in which R{x) is real. Specialization to real hermitian cases is not only useful for applications 
but can often make a more detailed analysis possible. 

A direct consequence of the fact that there is no first derivative term in ([1]) is conservation of the Wronskian W: 



A|^ = 0, W. 
dx 



(11) 



where '^-'^^w(x) and '^'^^u{x) are two exact solutions of Eq. Eq. pT|) also holds for complex values of a;, i?(a;), '''^^u{x) 
and (2)^(2;). 

A possible generalization of the Wronskian W to the vector case is 

Wn = (i)u(a;)-(2)u'(a;) - (2)u(x) • (i)u'(x). (12) 

The so defined Wronskian is conserved, if '^^'u(a;) and (2)u(a;) are solutions of Eq. ^ and the matrix R(x) is symmetric, 
again for complex x, R(x), (^)u(a::) and (2)u(a;). 
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For our purposes, another conservation rule will be useful, i.e., conservation of the generalized Wronskian Wn 
defined as follows, which holds if the matrix R(x) is hermitian: 



-^Wn = 0, WN = Re 
ax 



(«u(x),(2)u'(x))-((2)u(x),Wu'(:r))l, (13) 



where (^)u(a:) and ^^^u(x) are (complex or real) exact solutions of Eq. ([3]). 

If (i)u(x) and (2)u(a:) are two real vector functions, the two above definitions of the Wronskian are identical, and 
the current aj\j associated with the complex function 

u{x) = ^^'^u{x)+i^^^u{x), (14) 

is equal to the Wronskian Wn (= Wn)' 

aN = Im((i)u(x) + I (2)u(a;), ^^^u'ix) + i (2)u'(a;)) ^ Wn- (15) 

If these two vector functions are exact solutions of Eq. ([3]) for a real hermitian case, u(a;) is also an exact solution, 
and the conserved quantities ctat and Wn are the same. In other words, the current associated with any complex 
solution u{x) for a real hermitian case, is equal to the Wronskian Wn for '^^^u(a;) — Re u(a;) and (^)u(x) = Im u{x). 

In the scalar theory of PIA, it is often convenient or even necessary to go to the complex x plane. That is the case 
especially in higher orders, where the phase integrals are divergent at zeros and certain poles of R{x). If any such 
point is located on the real axis, it must be encircled in the complex plane. Going to the complex plane, however, 
in general violates the hermicity condition ([5]) or (fTU|) . Therefore in the hermitian theory of PIA, the independent 
variable will be assumed to be real. To remind the reader of this requirement, the independent variable in this paper 
is denoted by x rather than z. 



II. PHASE INTEGRAL APPROXIMATION IN THE SCALAR CASE 

The Phase Integral Approximation was introduced in 1966 by N. Froman'* in an attempt to improve the well known 
JWKB approximation. Later it was extended by introducing an arbitrary (base) function which could make the 
approximation work at its critical pointsi^ii^ Here, this approximation in its most general form will be rederived so 
as to make straightforward its generalization to vector cases. In this section, x and R{x) in Eq. ([1]) are allowed to be 
complex. 

In our derivation use will be made of a simultaneous transformation of the independent variable, a; — > xi, and the 
dependent one, u ui (Schwartz transformation), under which Eq. ([T]) conserves its reduced form: 



u'({xi) + Ri{xi) ui{xi) — 0. 



(16) 



The transformation x —f xi introduces the first derivative term in Eq. ([T|), which can be eliminated by using Eqs. ^ 
and ([6]) specialized to iV = 1. The result is 



ui = ql^"^ {x) u, qi{x) 



dxi 
dx 



Ri{x,)^q^\x){R{x) + S,[qi]}, 



(17) 



where is the nonlinear differential operator given by 



Sx[q] 



■,1/2 



-1/2 _ 3 
— 4 



q'{x) 



qix) 



2 q{x) 



_5_ 

16 



q'^ix) 



i<i 

4 q^{x) 



(18) 



[S-xVli] — , where {xi]x) is the Schwartzian derivative;-) Note that Sx[q\ is single valued \iq^{x) is, i.e., q{x) 

either single valued or defined up to its sign. 

Three properties of the operator Sx will be used in what follows: 
Homogeneity: 



5*0; [ckg] ~ Sx[q] if oi = const. 



Two rules for differencing the product: 



Sx[qi q2l] = Sx[qi] + 1 ^1 ^21 (^) ^ ^^[^^^j ^ ^^[^^j ^ g?5.j92l] 



xi — qi{x) dx. 



(19) 



(20) 
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and the third^ can be arrived at by considering three 

X2 {R '^^i?2) generated by 



The first two properties follows immediately from Eq. 
Schwartz transformations: 

X xi, as described above, xi — > X2 {Ri '^^-R2) generated by 921 = and x 

92 = ^ = 921 171, and requiring that (i)i?2 = ^^^i?2- 

Eq. can be used as an appropriate analytically solvable model for Eq. in which Ri{xi) should reflect basic 
properties of R{x) in some interval of x under interest. The problem then is to solve Eq. (|17p for qi{x), for given R{x) 
and Ri{xi). And with the proper choice of model, one can expect qi{x) to be smooth and "slowly varying" , and look 
for convenient approximation schemes. See Refs. 2, 9 and 11 for the lowest order theory, Ref. 14 for its generalization 
to higher orders and Refs. 15 and 18 for typical applications. 

If one is interested in solving the wave equation ([1]) in an adiabatic region, where the function R{x) changes very 
little on the characteristic scale of the solution u{x), the best model seems to be i?i(a:i) = c = const, e.g., c = 1. 
With this choice, ui{xi) — exp{±ixi), and Eq. (fT7|) leads to (we write q instead of qi): 



u{x) 



--(x) 



= q ^^"^{x) exp 



i q{x) dx 



(21) 



(22) 



where q{x) = [(7^(0:)]^/^ is double valued (defined up to its sign), q^{x) must satisfy 

q\x)-SM = R{x). 
and are constants. 

If q^{x) is an approximate solution of Eq. (j22p . the two functions u^{x) are also only approximate solutions of 
Eq. (IT|), called phase integral approximations. The name evidently appeals to the case of x and q{x) being both real. 
Another important special case is that of real x and q^{x) < 0, i.e., q{x) pure imaginary. In these two cases it is 
convenient to choose the constants and the sign of q{x) so that 



U^{x) = \q{x)\-'/^ < 



exp 
exp 



±i J \q{x)\dx'^ ifq'^ix) > 0, 
± / \q{x)\dx] ifq'^ix) < 0. 



(23) 



One should realize that the approximations pT|) and (|23|) (better or worse) always behave as if they were exact 
solutions of Eq. ([1]), e.g., they exactly conserve the current ai given by Eq. Q, and the Wronskian (fTTj) . Thus 
(cTi = Im[u±*(a:)u±'(a;)]) 



if q^{x) > 0, 
if q^{x) < 0, 



W = 



u'^ (x) u (.x) 
u+ {x) [x) 



-2i ifq'^{x)>0, 
-2 iiq'^{x)<0, 



(24) 



and W — zti 2c~^c~ in the general case of u{x) given by Eq. pT|) . 

If q'^{x) < 0, current conservation is actually a trivial consequence of the fact that u^{x) are real functions. A non 
trivial statement, however, is that the current associated with u~^(x) + iu^{x) is equal to —2, i.e., equal to W given 
by Eq. in accord with Eq. for iV = 1. 

All these nice features of (x) are due to the fact that these approximations are exact solutions of some equation 
of the form ([T]), in which R{x) for given q'^{x) is defined by Eq. (|22|) . This R{x) is single valued if q^{x) is, and 
is regular if q'^{x) is regular and non zero. Furthermore, R{x) is real for real x if q'^{x) is, which implies current 
conservation ([M)) . These facts were first pointed out in Ref. 16. And as both functions (PT|) or are solutions of 
the same equation ([T]), each linear combination of these functions (with real or complex coefficients) will also be a 
solution. Therefore it will also conserve the Wronskian W and the current ai . This property of two exact solutions 
is by no means obvious for two approximate solutions, in view of non-linearity of the Wronskian and the current. 

In Sec. IIVI the functions u^{x) will be generalized for the vector case to u*(a;) so as to conserve the generalized 
current Q in each approximation order. However, nothing analogous to Eq. ([T|) with R{x) given by Eq. (P^ . satisfied 
by u^(x), will exist there. Nevertheless, linear combinations of u^(a;) will be shown to also conserve the generalized 
current ^ in successive approximation orders. 

To construct a systematic approximation scheme for q{x) satisfying Eq. (|22p we assume that R(x) contains a small 
parameter A: 



R{x) ^ X-^G{x) +a{x), 0<A<1, 



(25) 



where G{x) ("greater" term) represents the dominant contribution to R{x) in the A limit and a{x) is an auxiliary 
function which can be chosen in any convenient way. In the scalar case and sometimes also in the vector case, the 
small parameter A can be eliminated from the final results by putting A = 1, see Sec. IVIII for more details. 



5 



Condition ((25)) is a quantitative statement expressing the adiabaticity of R{x). Indeed, if we freeze R{x) at its value 
for some x — Xq, the solutions of ([1]) will be 



u{x) 



exp{±i2Trx / L) if R{xo) > 0, 
exp(±27ra;/L) ifi?(a;o)<0, 



(26) 



where L — A27r|G'(a::o) + X'^a{xo)\~^^'^ is the characteristic scale for u{x) (the wavelength, or 27r times the e-folding 
distance). This scale is small as compared to that for R{x) (which is A independent). Hence, L can be expected to 
also be a characteristic local scale for the exact solutions. 

It is convenient to introduce A to Eq. (|2ip by replacing q ^ X^^ q. Finally we obtain, in view of 5*2, [A^^ q] = Sx[q], 
and with an appropriate choice of , 



u = q 



-1/2 



(x) exp 



iX ^ q{x) dx 



(27) 



Gix) - q^{x) + A2{5,[g] + a{x)} = 0. (28) 
Denoting q{x) in lowest order by Q{x) we obtain 

Q\x) = G(x). (29) 

This defines two solutions differing in sign, -iiQ{x), which can be improved in higher orders by adding terms propor- 
tional to A™, m = 1,2,..., 

(7(a;)== ^j/„(x)A'", yo{x)^±Q{x). (30) 

m=0 

Alternatively, we can multiply ±(5(a;) by one plus higher order terms: 

q{x) = ±Q{x)Y{x), Y{x) = ^™(^) ^o(^) = 1- (31) 

The functions Ym{x) are more convenient to deal with than ym{x), due to Yq{x) = 1 (in contrast to yQ{x) = 
■iiQ(x) ^ const). First, the applicability condition for the approximation in question if expressed in terms of Ym{x) 
is simply 

A'"|y™(a;)| « 1, m = l,2,.... (32) 

And secondly, the recurrence relations for Ym{x) are simpler than those for ym{x). An essential point is that equation 
for Y{x) which follows from Eq. ([^5]) is not more and in fact even less complicated than ([^5]) if the independent 
variable is appropriately changed. Using Eq. (|20[) we obtain 



Sx[{±Q)Y]^Sx[Q]+Q\x)S^[Y], (33) 

where 



C = ± J Q{x)dx. (34) 



Note that while x in Eq. ([T]) is usually a dimensional quantity (e.g., the space or time variable), C defined by Eq. 
is dimensionless. Finally, equation for Y{x) can be written 

(1 - Y'')Y' + X'{eoix)Y' + | [Y'{C)]'- ^^^"(0} ^ (35) 
where the additional term Sx [Q] coming from Eq. (j33p has been incorporated into another dimensionless quantity 

, ^ Sx\Q\ +a(a;) 
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The corresponding equation for q(x) is 

(g2 - q^)q^ + X^{a{x)q^ + | [q'{x)] ' - ^qq" (x)] - 0. (37) 

As is the only power of A occurring in Eqs. and ([57)l . Y and q can be expanded in powers of A^ rather than 
A. Thus, replacing in the expansions ((3T|l or ((30|) m ^ 2n, inserting them into Eq. |35|) or ([37)1 and equating to zero 
the coefficients of A^", we obtain the recurrence relations for or j/2n- Those for take the simple form {n > 1): 

^ r2„r2^- ^ Y2o.Y2pY2^Y2S + J] hi^2ar2/3 + fr^'JO^^S/S (C) " l>^2a>^2'MC)] = «• (38) 

n a+/3+7+(5— n a+/3— — 1 

Starting with yo(^) = 1, one obtains: 

>^2(x) = ieo, Y^{x)^~l[el + e'i{C)],---- (39) 
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Eqs. ((38)) and p9|) were first derived by N. Froman 
An explicit form of Yzn defined by Eq. (I38p is 



+ E hi^2ar2/3 + |l'2'a(C)5"2/3(C)- ^5^2.^2'^ (C)]], (40) 

a+/3=n-l 

where the tilde associated with the first two sums means that none of the subscripts a, (3, 7, 5 in these sums can reach 
the maximum value n, i.e., < a, /3, 7, (S, cr < n — 1. 

All functions F2„(x) are polynomials in p = 0, 1, 2, ... , with rational coefhcients, and the relevant formulas 
through 120(2;) were first obtained by J. Campbell^ by computer. 

General properties of Y2n{x) for arbitrary n are discussed in detail in Refs. 13 and 16. Note that all functions 12^(3;) 
can be expressed in terms of single valued quantities, Q^{x), eo(a;) and derivatives ^ of these functions. Furthermore, 
the relevant formulas contain no complex coefficients, see the definition of eo(a;), Eq. (|36[) . and the identities 

yLiOy^^O = Q-\x)YiJx)Yip{x), Yi'piO = Q-%x)[Yi'p{x) -\Q-\x)Q^'{x)Y^p{x)\. (41) 

As a consequence, all functions Y2n{x) are invariant under change of sign of Q{x) and are real if x, Q^{x) and a{x) 
are real {Q^{x) > or Q'^{x) < 0). Taking m — 2n in the expansion ([3T|) and truncating it at n = A/", we obtain 

^^ 

q{x) = q2M+i{x) = ±Q{x) ^ Y2n{x) A^". (42) 

Inserting this q{x) into Eq. (|27|) . we obtain two linearly independent approximate solutions of Eq. ([T]). They are called 
phase integral approximations of order 2M + 1 (as they are related to the JWKB approximations of this order) . We 
recall that x and Q'^(x) can be complex, but specialization to real values is useful for applications, see Eq. (|23p . 

The recurrence relations for yn{x) equivalent to (|38p can be obtained from Eq. (|37p . if a(x) is expressed in terms of 
£0(2;) and Q^{x) by using Eqs. and The result is 



+ H { (^'^0 " Sx[Q])y2o.y2p + ^y2J^)y2pi^) - \v2c.y%{x)\ = 0, (43) 

which is evidently more complicated than Eq. ([55)1 . Additional terms with Q'{x) and Q"{x) present in are 
necessary to cancel similar terms produced in y2n{x) when differencing Q^(x) present in Eq. (|43p . Only after these 
cancellations, the actual dependence of y2n{x) on eo(C) can emerge (rather than separately on Q'^{x) and a(a;)). This 
type of dependence is directly seen from Eqs. ([55)1 and ([55)) - (|in|) but is rather hard to see when starting with Eq. (|37p . 
Note also that y2n{x) (= Q{x) 1^2n(a;)) is complex if Q^{x) < 0, in contrast to Y2n{x) which is real if Q^{x) is real. 

A relatively simple program in Mathematical'^ based on Eq. (UOl) enables one to generate the corrections Y2n{x), 
n = 1,2, . . ., either in their general form analogous to Eq. ([55)1 (with possible transformation of the derivatives ^ 

into ^), or for each given choice of R{x) and a{x). 
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III. GENERAL THEORY OF THE PIA IN VECTOR CASE 



In the scalar case, the starting point of the phase integral theory are Eqs. (p5)) and ([27]) • Their generalization to 
vector cases is straightforward, i.e., 

R(x) = A-2G(a;) +a(x)I, < A < 1, 



u(x) = s{x) q ^^^(x) exp iX ^ q{x) dx 



(44) 
(45) 



where I is the unit matrix, s(x) G Ti,^ and a{x) is an auxiliary function. If we choose a{x) = 0, Eqs. and ([IS]) 
become equivalent to those proposed by Fulling.^ 

Inserting Eqs. (|44)) and ([45]) into Eq. ([3]) and multiplying by A^, we easily find an equation that governs the new 
unknowns s{x) and q{x): 



(x) - q^{x)l\ •s(x) + X2iq{x)s'{x) + \ s" (x) - s'(x)^^ + f S,[q] + a{x))s{x) 
i I q[x) V / 



0. 



(46) 



In Fulling's theory, the matrix R(a;) was assumed to be hermitian. This is equivalent to the hermicity of the matrix 
G(a;) now entering Eq. if we assume that a{x) is real for real x. 

Note that the number of new unknowns (A^+ 1) is greater than the number of Eqs. (|46|l (N). Therefore a constraint 
upon the unknown vector s(a;) is needed to guarantee the uniqueness of u{x). 

All theories developed in this paper will start with Eqs. (I44|) -(l46 |) . They will differ in the adopted form of the 
constraint. 

For = 1, on replacing s(x) 1, G{x) G{x) and 1^1, Eq. reduces to and s{x) is A independent. 
Therefore in the general case of > 1, the expansion of s{x) in powers of A must start with the A independent term: 

S{x) = ^ra{x)X^. (47) 

Using Eq. pS)) in lowest order (A°) and denoting, as in the scalar case, the q{x) in lowest order by Q(x), we obtain 

G(.T)-g2(a;)l].so(a;) = 0. (48) 

This indicates that Q'^{x) must be an eigenvalue of the matrix G(x), and So(a;) is the corresponding eigenvector. If 
G{x) is hermitian, the eigenvalue Q^{x) is real i.e., Q{x) is either real (if Q^{x) > 0) or pure imaginary (if Q^{x) < 0). 
Fulling's paper— was restricted to Q^(x) > 0. Here, as in the scalar case, both situations will be discussed. 
The eigenvalue Q'^{x) can be found as the solution of the characteristic equation 



det 



G(x) -g'(x)i 



= 0. 



(49) 



The LHS of Eq. (gHl) is a polynomial in of degree N. 

For reasons explained in Sec. [TTl it is convenient to assume that the expansion of q{x) in powers of A has the 
form (pij) . Repeating the arguments following Eq. (PT]) . where the 0th order equation (I29p must now be replaced by 
Eq. (HS)) . and expressing the derivatives Qj^) in Eq- 



in terms of 4? we obtain 

dC, 



Y^{Q-^G - YH) ■ (s - So) + (1 - Y^)Y^so + Xi2Y'^s'{0 



+ x^[y\"{0 - ry'(C)s'(C) + [^^y^ + l[y\C)f - \YY"{o] s} = 0, (50) 



where C, and Cq are defined by Eqs. ([M)) and ([55]) . A distinguishing feature of this equation as compared to that in 
the scalar case ([55)1 (obtained if we replace s, So ^ 1) is the presence of the A term. It contains the pure imaginary 
coefficient i and changes sign if Q is changed into — Q (i.e., dC, —dQ- This term is responsible for differences 
between the vector and scalar theory. 

The LHS of Eq. (f50|) is a polynomial in A containing only positive powers A™, m — 1,2,... (as s = sq and Y = 1 
in lowest order). Equating to zero the coefficients of A™ we obtain: 



Q- 



E 



Y^YpY^Yss, 



a+/3+o"— m 
cr>l 



cr>l 



a+/3— m 



Y Y^YpY^Ys) So + ^ -0, 

a+/3+7+i5— m m — 1,2 



(51) 



8 



E 



2i 



E 

Q+/5+74-cr— m — 1 



(52) 



where the second sum on the RHS in Eq. ()52p must be dropped if m = 1. Writing down exphcitly terms in the 
first four sums in Eq. (|5ip corresponding to maximal values of a, (3, 7, 5, a (= m), this equation can be written as an 
implicit form of the recurrence relations for Ym and Sm {m > 1): 



F™so-iQ-"(G-Q^I)-s,„-b„ 



(53) 



where b,„ depends on and Sg- with a, a < m — 1. Using Eqs. (|5T |) ~ ([55)) . we can easily find the recurrence relations 
for b,„ (m > 2): 



2[ E 

a+/3+cr— m 
cr>l 



y„y^ (s, + 2(i;so - b,)) - ^ Y^YpY^Yss, 



CT>1 



51 Y^YfjY^Ysjso+ J2 

a-\-j3—m o+/3+7+(5— m rn— 1,2 



(54) 



where the tilde associated with the first four sums means that none of the subscripts a, (3,j, S, a in these sums can 
reach the maximum value m, i.e., < a, /3, 7, (5, cr < m — 1. The process of generating b2, bs, . . . from Eq. (|54p can be 
programmed in Mathematica;^ and bi can be found by dropping the sums with tildes and with a + P + a = m— 2. 
The results are: 



bi = is'oiO, 

b2 = i(eo-n2)so + iFis;,(c) + is;,'(c)-Fisi + is;(c), 

b3 - -[rly2 + ir/'(c)]so + [^r2-ir/(c)]s^,(c)-i(rl2 + 2y2-eo)sl 

+iYis[{C) + is'/(C) - Y1S2 + z s[,(C), . . . , 



(55) 



= hm YiS 



m — 1 



IS, 



rn — l 



(C), 



(m > 2), where hm depends on Yq and So-, sJy(C), . . . , with a < m ~ 1 and cr < rn — 2 (hm is independent of Sm-i, 
s;,_i(C),...). 

In the original treatment)^ Fulling expands p{x) (= 9^(2:)) rather than q{x) in powers of A (= u ^ in his notation, 
u > 1): 



Pm{x)X'' 



(56) 



This nonstandard approach introduces unnecessary complication, due to the appearance of ^yp{x) (^— q{x)^ in the 
integrand in Eq. (|45| and in the A term in Eq. (|46| . As a consequence, evaluation of the phase integral (|45p becomes 
nontrivial even in the simplest cases, and it would be much more difficult to find equations analogous to (|53p and 
(|54p with Pa instead of Yq, valid for any to. The functions p,„(x), if needed, can easily be expressed in terms of 10(2;) 
(= l),Yi{x), ...,Ymix): 



Pm{x) = Q'^{x)YYa{x)Ym-a{x), TO > 0. 



(57) 



In the general discussion that follows, the assumption of x and a{x) being real and R(a::) hermitian is not necessary, 
unless a special comment is made. 

An important role in vector theory is played by single valued quantities (invariant under the replacement Q — *■ — Q, 
i.e., d( —d(): Q'^{x), eo{x), G{x), So{x) and derivatives ^ of these functions. Using Eq. ([53]) along with (l55|) it 
can be seen that bi, Yi and Si are double valued, b2, Y2 and S2 are single valued etc. In general, the even order 
corrections remain invariant under the replacement Q — > —Q, and the odd order corrections change sign. This means 
that, in analogy to the scalar case, the even order corrections can be expressed in terms of single valued functions 
Q^{x), eo{x) etc. Using Eqs. ([3T]) . (|45| and ((47|l . we obtain two phase integral approximations u^{x): 



u{x) — u^(x) = c'^s^lx) q''^{x) 



exp 



i X / q (x) dx 



(58) 
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where c are constants and 



(59) 



Note that the even order contributions to q^{x) are double valued (like those in the scalar case), whereas the odd 
order contributions (specific to vector case) are single valued. 

Eqs. ([58|) ~([59 |) generalize ((2T|l and ([3TI1 (with m — 2n) to the vector case. The main difference in comparison to the 
scalar case is that u='= (x) are not solutions of the same differential equation of the form ([3]) . 

In hermitian vector cases, where Q'^{x) is real, we can take 



±Q{x) 



±\Q{x)\ ifg2(a;)>0. 



Ti\Qix)\ ifg2(x)<0, 
and choose the constants so that (^q^{x) — \Q{x)\Y^{x)') 



(60) 



u^{x)=s^{x) q^{x) 



-1/2 exp 



exp 



±iX-^ J q'^ix) da;J if Q^{x) > 0, 
±A^i / q^{x)dx] iiQ^{x) < 0. 



(61) 



Note that, in general, So(a;) is complex. Definite statements concerning reality etc. can only be made in real hermitian 
cases in which Sq{x) is also real. 

Thus in real hermitian cases, if Q^{x) > (as in the original Fulling's theorj*®), it can easily be seen by inspection 
that the even order quantities, b2„, and S2„, n — 1,2,..., are real, and the odd order ones, b2„_i, etc., are 
pure imaginary, see Eqs. (|53l) - (l55|) in which d( = ±|g(a;)| dx is real. And if Q'^{x) < 0, d( = =F« |g(2;)| dx becomes 
pure imaginary which makes both the even and odd order corrections real. In any case, S2n+i(a^) and y2n+i(a;) (pure 
imaginary if Q^{x) > and real if Q^{x) < 0) correspond to the upper sign in Eq. ((60)) . Note that in both cases, 
the single valued contribution to the integral coming from Y2n+i{x) is real. This contribution slightly modifies the 
amplitude of the phase integral approximation in higher orders. Its role is similar to that of the factors in the first 
line in Eq. (j61[) . The main x dependence of the approximations u^{x) is given by the contribution to the integral 
coming from Y2n{x). The corresponding exponential (like that in the scalar case) either exhibits strongly oscillatory 
behavior if Q^{x) > 0, or exponential growth or decay if Q^{x) < 0. Note that if Q^{x) > 0, we obtain in analogy to 
the scalar case: 



-{x)=[u+{x)Y 



(62) 



And if g^(x) < 0, the approximations u^(a:) are real functions, again in analogy to the scalar case. 

When solving Eq. ([53]) for Ym and s™, an important point is whether Q^{x) is a degenerate eigenvalue of G 
{d > 1) or non degenerate {d — 1), where d is the dimensionality of the linear subspace of eigenvectors so{x) 
corresponding to the eigenvalue Q^{x). Denoting by Ti.^ the N — d dimensional orthogonal complement of TC^ and 
introducing orthonormal bases: in H'^, {cfc}, k — 1,2 ... ,d, and in H"*", {e^}, fc = 1, 2, . . . , — d, we obtain 



{ej,ek) = Sjk, E(efe,s)efe = s if seH'^, 



(63) 



k=l 



N-d 

{el,ei) = 6,k, 5](e^s)e^^s if seH^. (64) 

k=l 

Each vector s e can be decomposed into its component belonging to H'^, Ps, and that belonging to the orthogonal 
complement, s^'^: 



s = Ps + s^, Ps^=0, 

where P is the orthogonal projection operator, acting in Ti^ and projecting onto W^: 

d 

Ps = Y{ej,s)ej. 



(65) 



(66) 
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Using Eq. ([65]) for s = s„j in Eq. (|53|) . the contribution to the LHS coming from Ps,„ will be zero, i.e., we can replace 
in Eq. ([55)1 s„ — + s^. Multiplying the result by e^, we obtain equations governing the coordinates of s:^, (e^,s^). 



in the basis {e^}: 



N-d 



^(G|,-Q%fe)(e^,s^J--2Q2(e|,b„), Gf, ^ {ef ,G ■ ei), (67) 

k=l 

j, k — 1,2, . . . , N — d. The linear set of — d equations (j67|) is nonsingular, and its solution defines the orthogonal 
component s^: 

N-d N-^d 

si = -2Q2 ^ 4fe(e^b„)> ^ (G^ - Q^l^r\ (68) 

j=i fe=i 

where G"*- is given Eq. ((67)) and is the unit vector in H.^ . 

Equations involving coordinates of Psm are obtained on multiplying Eq. I|53p by the basis vectors ek ■ It is convenient 
to choose the basis {ek} so that Sg is directed along one of the basis vectors, e.g., 

So = |so| ei. (69) 

Multiplying Eq. (j53p by Sq we obtain an explicit form of the recurrence relation for Y„i: 

y,n - |sor'[i |so| (ei, G • s^) + (so,b„0], (70) 

and the remaining products lead to 

\Q~^{Gk,G-si) + {ek,h„,) = Q, fc > 1. (71) 
Using Eqs. ([M|) - ([7T|) for m = 1 and recalling ([55)1 we obtain 

Yi = |sor2[iQ-2|so|(ei,G.s^) + zQ-i(so,s;,(a;))], (72) 

\Q-^{ek,G-si) + iQ-^{ek,s'^{x)) =Q, fc > 1. (73) 

The following discussion in the next two sections will be given separately for the hermitian and non-hermitian G 
matrices. 

IV. HERMITIAN THEORIES OF THE PIA 

In this section we assume that the G(a;) matrix is hermitian which in general requires x to be real. This assumption 
simplifies the theory as in that case 

{ek,G-si)^0, fc = l,2,...,d, (74) 

(G • efe — ek). Furthermore, the matrix G-'- defined by Eq. (p7|) is hermitian. 
Eqs. dill) and ^ take the form 

Y,^tQ-'\so\-^So,s'o{x)), (75) 

(efc,s;,(a;)) =0, fc > 1. (76) 

It can be seen that Sq{x) must be orthogonal to all basis vectors which are orthogonal to Sq. In view of these d — I 
requirements one might be tempted to assume that Sq(x) is also orthogonal to ei, 

(so,s;,(a;)) = 0, (77) 

i.e., Sq(x) e H^. That is because in that case, Yi given by Eq. ([75)1 is zero, 

Yi{x) = 0, (78) 
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which means an essential simphfication of the theory, see Eqs. (|55|) . And as Eq. ((77)) impHes (so,So) = const, then 
with this assumption, the simplest choice is 

(so(x),So(a;)) = 1, i.e., Sq = ei. (79) 

The constraint 

(ei,e'i(a;)) = 0, (80) 

is often fulfilled automatically by an orthonormal basis. Otherwise, it can always be fulfilled by an appropriate choice 
of the phase factor in ei,® i.e., if we take 



ei = exp(?; 6?i) Gi , 9i = i J {ei,e\{x)) dx, 



(81) 



where 6i is a real quantity, and {e^} is an arbitrary orthonormal basis in "hC^ (e^ = e^, if fc > 1). 
Eqs. dZHl) and ^ lead to 

pi=0, p2^2Q^ix)Y2{x), p3^2Q^{x)Y3{x),... (82) 

Eq. ([70]) for Y^ becomes 

r„ = (ei,b„), (83) 

where is given by Eq. ([55)1 simplified by Yi{x) = 0. 

Note that the RHSs of Eq. (|53|) with bi, b2, etc. given by Eq. ((55)) with Yi{x) = 0, are simpler than the analogous 
expressions given by Fulling, see the Appendix in Ref. 8. That is because we are expanding Y rather than q in powers 
of A, where Yq = 1 (in contrast to yo = Q 7^ const) and where the C variable is the natural choice. Working with this 
variable is very convenient for derivations and presentation of final results. 

Replacing in Eqs. ((69)) - ((74)) m ^ m + 1 and using the last of Eqs. ((55)) along with ((78)) and ((79)) . we obtain 



Y„,+i^{ei,h^+i)+iQ ^ {ei,s'^{x)), (84) 



{ek,s'^{x)) = iQ {ek,hm+i), fc > 1, (85) 

where b^+i depends on Yq, with a < m and on So-, s'^{x), etc. but with a < m ~ 1. This means that Eq. ((85)) is the 
actual and only constraint upon s„j that follows from Eq. ((53)) for given m {— 1,2, . . .). Using the decomposition ((65)) 
we can write 

{ek,s'„^{x)) = (^Gfe, -^s,-|,J + (efe,-^(Ps,„)), (86) 



for any fc = 1, 2, . . . , d. In order m, the first term on the RHS is known, see Eq. ()68p . and in the second term, the 
derivative ^ can be shifted in front of the scalar product if the basis is appropriately chosen. Indeed, using the 
definition of the projection operator P, Eq. (|66p . we can write: 

d d 

— (efc,Ps„) = (ek,—(Ps,n)^ +^(ej,s„)(efc(x),ej), (87) 

where each term in the sum over j in Eq. ((57)) is zero if {efc(x)} is the Kato basisr' characterized by (j, fc = 1, 2 . . . , d) 

ie,ix),e'k{x))^0, i.e., Pe'kix)^0. 



Hence, in the Kato basis (which in particular satisfies our earlier requirements: (|76p with Sq = ei and ()80p ) we can 
determine the coordinates (efe,Ps,„), fc > 1, by using Eq. ((55)) and ((57)) along with ((55)) in Eq. ((55)) and integrating. 
The result is 

(efe,Ps„) = / (efc, i Qb,„+i - s^/(a;) ) dx, fc > 1. (89) 
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These coordinates are thus defined uniquely by the compatibiUty condition in order (m+1), Eq. ([85)) . At the same time 
the coordinate (ei,PSm) measured along ei is an unspecified function. The choice of this function leaves unchanged 
the remaining unknowns in order m, i.e., and the coordinates (e^, Ps™) measured along e^, fc = 2, 3 . . .. However, 
this function will influence Ym in next order, given by Eq. ([M]) . Using again Eqs. ([55]) and ([57| in the Kato basis now 
for fc = 1, we obtain 

Ym+i = (ei, hrn+i+iQ~^-^si^ + i ^ (gi , Ps™) . (90) 

The choice of the coordinate (ei,Ps,„) must be compatible with general properties of the even and odd order 
corrections discussed in the previous section. This will be the case if we take 

(ei,S2„) = f2n{Q'^{x),eo{x), ...), (ei,S2„-i) = iQ{x) /2„-i (Q^(x), eo(x), . . .), (91) 

n = 1,2, . . ., for any functions f2n and /2n-i of the mentioned earlier single valued quantities and their derivatives. 
These functions must be real for real arguments. 
Note that Eq. feS]) implies 



Re{e,{x),e',{x))^0. (92) 

Hence if lm{ej{x),e'i^{x)) = 0, as is the case for real eigenvectors, each orthonormal basis is the Kato basis. 
Also in complex hermitian but non degenerate vector cases (ImG ^ 0, d = 1) the orthonormal basis ([63)1 satisfying 
Eq. ([80|l is the Kato basis. Only in complex hermitian and degenerate vector cases (ImG 7^ 0, d > 1) must one solve 
the nonlinear set of 1st order ODEs for the coordinates eji {x) following from Eq. ([55)1 

N 

J2e*i{x)e',iix)^0, j,k = l,2...,d. (93) 
/=i 

The results of this section pertaining to hermitian vector cases can be summarized as follows. Given the orthonormal 
bases: {e^} in and {e^} in Ti.-^, where {efc} must be the Kato basis satisfying Eq. ([SS]) . an explicit form of the 
recurrence relation ([55)) is given by Eqs. ([S5|) . ([59)) and ([55)). which define s;^, Psm and Ym, respectively, in terms 
of the same quantities in lower orders. The only freedom in each order is the choice of the coordinate (ei,Psm). 
A characteristic feature of the vector theory, in contrast to the scalar one, is the presence of integrals in recurrence 
relations. However, they are only present in degenerate cases {d > 1). 

We can try to choose the unspecified coordinate (ei,Psm) = (ei,Sm) in each order so as to get the PIA in the 
vector case as close as possible to that in the scalar case. We recall that in the hermitian vector cases, exact solutions 
of Eq. ([3]) conserve both the generalized current aN and Wronskian Wn given by Eq. ([T3)) . In the scalar case, both 
these quantities are conserved also by the PIA. In vector cases, the freedom in choice of the coordinate (ei,Sm) can 
be used to conserve one of these quantities. The second one will be seen not to be conserved. Nevertheless, in real 
hermitian cases both the current and the Wronskian can be conserved in view of Eq. (|15p . 

If Q^{x) > 0, it is better to conserve the current ctat for each of the complex PIA u^(a;) rather than Wn for these 
two vector functions. That is because if this choice is made in a real hermitian case, two real functions Reu*(j:) 
and Imu^(x) (which are approximate solutions of Eq. ^) will conserve their Wronskians in view of Eq. ([T5)) . And 
if Q^{x) < 0, it is better to conserve Wn for u+(x) and u~{x) rather than the currents ajy for each of these 
approximations. These currents will not in general be conserved (for complex u^{x)). However, in a real hermitian 
case (where the real approximate solutions u^(a;) conserve the current identically), the current ctat associated with 
an approximate complex solution u+ (x) + i u+ (x) will also be conserved, again in view of Eq. ([T5)) . The first scenario 
was put into practice by Fulling^ and the second one will be be described in this paper. 

In the scalar case, conservation of the current af and Wronskian W (both involving the products of u^{x) and 
u^'{x)) is due to the fact that the x dependent factor \q{x)\ coming from u^'{x) is multiplied by two factors |(7(x)|^^/'^. 
This mechanism will work also in hermitian vector cases, see Eq. ()6ip . if Y^{x) is real and positive. Thus, assuming 
that Y^{x) > in successive orders (which must be checked a posteriori) and leaving out the superscripts ± we 
obtain ((Tjv = Im (u, u'{x)) , q{x) = \Q{x)\ Y{x)) 

^ fg-i(x)Im(s,s'(x))±A-i(s,s) if Q^x) > 0, 

\q-^x)lm{s,s'{x)) eiiv[±2X-^ Jq{x)dx] ifQ\x)<0. ^ ' 

This quantity will be conserved in each order if two constraints proposed by Fulling^ are fulfilled in each order: 

(s, s) = l, normalization, (95) 
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(s,s'(a-)) =0. 



(96) 



In fact, the actual constraint is Eq. ((96l) . whereas ([95]) must be fulfiUed at some fixed value of x. That is because 
Eq. ([96]) implies (s,s) = const, and so this constraint is enough to guarantee cjn ~ const. The second constraint ([95]) 
is rather cosmetic. In this paper it will be fulfilled in lowest order only. 

In lowest order, the constraints ([95)1 and (|96p reduce to Eqs. ([7^ and ([50]) . In higher orders Eq. (|96p takes the form 



={), TO =1,2. 



(97) 



(which implies X]a=o(^a' ^ const). 

If Q'^{x) < and Y^{x) > 0, Wn for u^{x) contains the factor exp [X-^ J \Q{x)\{Y+ {x)-Y- (x)) dx]. Conservation 
of this factor requires that y+(a;) = Y^{x), i.e., 



r2„+i(a:)=0, n = 0,l,. 
With this condition fulfilled, Wn is given by {Y = Y+ = Y- > 0) 



Wn= [\Q{x)\Y{x)] 'Re[(s+,s-'(a;))-(s-,s+'(x))j [(s+, s") + (s", s+) 



Eq. ([59]) implies 



(s±,sT(x)') = ^(TA)"^(-l)"(s„,s:„„Jx)). 



(98) 



(99) 



(100) 



Thus if 



^(-l)"(s„,s:„_Ja;)) = 0, m = 0,l,2,. 



(101) 



a=0 



we obtain 



(s+,s-'(a;)) = (s-,s+'(a;)) =0 

in successive orders. This in turn leads to 

[(s+(x),s-(a:)) + (s-(x),s+(x))]'-0. 

Thus, if Y2nix) > and conditions ([98]) and ([TOT|) are fulfilled, Wn is conserved in successive orders. In lowest order, 
Eq. ()10ip is satisfied in view of Eq. (sq = ei). 

If Q^{x) > and Y^{x) > 0, conservation of Wn for u^{x) also requires Y2n+i{x) to vanish. With this condition 
fulfilled, Wn takes the form 



Wn= [\Q{x)\Y{x)] ^Rc <^ (s+,s-'(a;)) exp -2A"i \Q{x)\Y{x)dx 



- s ,s 



-'{x)) exp 



2A- 



\Q{x)\Y{x) dx 



(102) 



This will be zero (and thus conserved) in successive orders, if condition (llOip is fulfilled. 

An important point is that either condition ()97p or (llOip with to > define uniquely the unspecified coordinate 
(ei,Sm) (up to the integration constant). The conditions in question can be written (5(to = 1) = 0): 



{e,,s'^{x)) = -{±ir {sm,e[{x)) - S{m), 5(to > 1) ^ ^ (±l)"(sX.-a(^))- 
Using this result in the identity 



(103) 



a=l 



(ei(x),s,„(x))' = (e;(a;),s„) + (ei,s^^(x)), 



(104) 
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and integrating we obtain (see also Eq 

(ei,s„)= / {e[ix),si) - {±ir{e[{x),siy ~ Sim) 



dx. 



(105) 



Grouping together in S(m = 2n) and S'(m = 2n + 1), the first and last term, second and last but one etc., and picking 
up total derivatives we end up with (n = 1, 2, . . .) 

(ei,S2„) = z2Imy'[(e;(x),s^„(x))-(±iri(s„(x),s'„(x)) 

n-l 

-;^(±l)"(s2„-„(2:),s^(a;))] dx - {±1)^ ^{sJx),Sr,{x)) 

a=l 
n-l 

-^(±l)"(s„(:E),S2„-a(x)), (106) 

a=l 

where the sums over a must be dropped li n — 1, and 

. n 

)= / /(x)dx-^(±l)«( Sa{x),S2n+l-a{x)), (107) 



/(^) 




ei(2;),sij+i) + ^(s^(a;),S2„+i_a 

a = l 
n 

ei(a;),si-„+i) + ^(-l)"(s^(x),S2„+i 



(108) 



The last two equations are also valid for 71 = if sums over a are dropped. In that case, the upper line in (|108p 
gives the FuUing's result^ if we choose a{x) = 0. In general, for any n > 0, the upper line in Eq. (|108p and upper 
signs in (|106p and (|107p refer to the current conserving theory and lower ones to the Wronskian conserving theory as 
developed in this paper. 

Note that in a real hermitian case, the integrand in Eq. (|106p is real both if Q^{x) > (scalar products of either 
real or pure imaginary factors) and if Q^(x) < (all factors real). Therefore in that case, (ei, S2„) contains no integral 
contribution in either current or Wronskian conserving theory. The integral contribution is absent also in (ei,S2„+i) 
in the current conserving theory if Q^{x) < and in the Wronskian conserving theory if Q^{x) > 0. 

In applications, the matrices R(a;) and G(a;) are often real and symmetric ("real hermitian case"). In that case, if 
Q^{x) > (as in the original FuUing's theory^), the odd order corrections, l2„_i and S2,i_i, n — 1, 2, ... , are pure 



imaginary, see the text following Eq. ((6T|) . In view of Eq. ([57)) the same is true of P2n~i- Confronting this fact with 
FuUing's statement^ (without proof), that in the hermitian vector case with Q^{x) > 0, all corrections Pm(a;) are real, 
we would obtain P2n-i = (i.e., l2n-i = 0) for real hermitian cases. We were unable to prove this fact in general, 
but the examples given in Sec. I Villi do confirm this behavior. Note also, that in real hermitian cases with Q'^{x) > 
and (ei, s™) given by Eqs. (jl06p - (|108p (upper line and signs), not only u^(a;) conserve the generalized current in each 
order, but that is the case also for linear combinations of u^(a;). Indeed, introducing w(a;) = c+u+(x) + c~u~(a;) 
and denoting by crjv, cr^ and cr^ the currents associated with w(x), u+(a;) and u^(a;) we obtain 



N 



\c-?a 



N 



Im 



(u+,u-'(x)) +c+c-*(u-,u+'(x)) 



(109) 



where Im[ ] = in view of Eq. (|62p . At the same time, the current associated with linear combinations of two PIAs 
generated by two eigenvalues (5^j^^(a;) and (3^2) (2^) in FuUing's theory in general is not conserved. For example, for 
zero order PIAs we obtain (if Q'^^^-^{x) ^ Qj^-^^x), which imphes (ei (1), ei (2)) = 0) 

(^N = |c(i)|VAr(i) + |c(2)|VAr(2) +2Im c^i)/(\)(x)c(2)/(2)(a:;)(ei(i),e'i(2)(a;)) , (HO) 
/(i,2)(a^) = |Q(i,2)(a;)r'/'exp[±z f \Q(^,^2){x)\ dx] , (111) 



where in general Im[ ] 7^ const. 
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Eqs. (|68|) . ((84| and ([89]) along with (|106p - (|108p are the recurrence relations in explicit form, from which the higher 
order corrections and can be determined for any m > 1. They confirm our earlier observation that in the real 
vector case (where can be assumed to be real) and Q'^{x) > 0, indeed are the even order corrections S2q and 120: 
real, and the odd order ones pure imaginary. It can be shown that Y2 (= ^ P2) is also real in complex hermitian 
cases, see Eq. (32) in Ref. 8. 

For the real hermitian cases with Q^{x) < 0, Q{x) becomes pure imaginary, but all other quantities in the recurrence 
relations, including and hm+i given by Eqs. (j55p and s;^ given by Eq. (j68p . are real. This implies (as mentioned 
earlier) reality of the functions u^(a;). 

Note that the case oi d ~ N (full degeneration) is trivial, as in that case each vector u S is an eigenvector of 
the G matrix corresponding to the eigenvalue in question. This means that G • u = Q'^{x) u, and Eqs. ^ and 
(|44)) lead to 

u"(x) +i?(a;)u = 0, (112) 

with R{x) = X^^ Q'^{x) + a{x) [N identical scalar cases). 

For rf = — 1 , the orthogonal complement H.^ is one dimensional and the orthonormal basis in H.^ contains only 
one vector e-'-. Denoting by s-^ any vector belonging to Ti^ we obtain = |s^|~^ s-*- and Eq. leads to 

= -2Q2^-i[(s^,b„0-i^i(s^,s,„_i)+zQ-i(s^,s:„_i(x))]s^, 
D = (s^,G-s^)-|s^|2q2, (113) 

(The first two terms in square brackets in pi3p must be dropped if m = 1.) Thus in that case, finding is very 
simple, and the main algebraic problem is to determine the Kato basis in T-C^. Furthermore, in each order one has to 
find d integrals, see Eqs. and (fT05ll -(fT08 | . 

The simplest alternative to calculating (ei,s„i) from Eqs. (|105p - (|108p is to choose f2n — f2n-i = in Eqs. ([9T|l . 
leading to 

(ei,s,„) = 0. (114) 

This simplifies the theory by eliminating an integration in each order. This choice, which will be referred to as 
"simplified hermitian theory" , seems to be especially attractive in the non-degenerate case in which the eliminated 
integration is the only one present. The price, no current or Wronskian conservation in higher orders, is probably 
worth paying in view of advantages of this theory, i.e., its simplicity and absence of logarithmic terms in the higher 
order corrections, see Sec. IVIIII 



NON-HERMITIAN THEORY OF NON-DEGENERATE VECTOR CASES 



If the G matrix is non- hermitian, the main complication is non vanishing of the product in Eq. (|7^ . As a 
consequence, Eq. (|73p cannot in general be fulfilled. That is why our non-hermitian theory in this section will only be 
given for the non-degenerate cases (rf = 1) in which Eq. (|73p is not present. In this section, the independent variable 
X can be complex. 

In general, Yi given by Eq. ([7^ will be non zero. This will introduce complication in Eqs. ([55)) . ([55]) and ([70)1 . In 
this connection there is no need for normalization of Sq or imposing the differential constraint (|77p . Choice of the 
multiplication factor (or function of x) in the eigenvector should make the theory as simple as possible. This should 
be the only criterion also in those hermitian cases in which the integral in Eq. (j8ip is not expressible in elementary 
functions, again leading to Yi{x) ^ 0. 

It will be assumed that the only component of Ps,„ which may be non- vanishing for d = 1, i.e., that measured 
along So, is zero: 



(so, Sm) 



N 
J = l 



0. 



(115) 



This means that Sm is equal to s^j, defined by Eq. ([68|) . and Ym is given by Eq. (|70p. both in terms of the — 1 basis 
vectors e^. Dealing with these basis vectors is inconvenient and can easily be avoided as follows. 

Due to our assumption of d = 1, the rank of the matrix G — Q^I in Eq. ([53)1 is — 1. Assuming that 



G22 — Q 
G32 

Gn2 



G23 

G33 - 



G 



N3 



G 



G2N 
G3N 

NN 



^0, 



(116) 
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and introducing (denoted by a bar) vectors and matrices in the — 1 dimensional Hilbert space Ti,^ ^ {j, k — 
l,2,...,iV-l) 

Srak=Srak+\: bmk = bmk+l, §{0:, 1) k = G k+1 1 , g{l, a) k — Gi k+1 , Gjk~Gj+ik+l, (H^) 

we can write Eqs. ([55)1 and (IllSp in vector form: 

(Gil - Q')s™i + g(l, a) • s„ = 2 Q2 (^Ym - 6ml), (118) 

(G - g^I) . s™ = -s™ig(a, 1) + 2 Q2 go - b„), (119) 

Sqi Sml + (so, s„,.) = 0. (120) 
Using Eqs. pisp multiplied by Sq^, (|119p multiplied by \soi\^ and (|120p . we can easily find 
Sm = 2(5^SoiD"^ • (6,„iSo - soibm), 

D = \soi\''{G-Q^i) + iGii-Q^)sos*-s*,Sog{l,a)-.soigia,l)s*. (121) 

The last three terms in D given by Eq. (|12ip contain dyadic products of vectors in H^^^, defined as (pq)jfe — PjQk- 
Now Sml can be determined in terms of by using Eq. (|120p . 

Sml = -(sO, Sm/Soi), (122) 

and Yjn can be found from Eq. (|118p . In the hmit sqi ^ 0, s,„/sqi in Eq. (I122p is well defined, see Eq. I|12ip . This 
will be illustrated by examples given in Sec. IVIIII 

VI. PHASE INTEGRAL APPROXIMATION FOR N = 2 

The case of two equations ([2]) {N = 2) is the only situation in which the algebraic part of either the current or 
the Wronskian conserving theory (with d = 1) in higher orders is very simple. This is because in that case, both the 
eigenspace Ti,'^ and its orthogonal complement H'^ are one-dimensional. In any of three hermitian theories and in the 
non-hermitian one developed in Sees. IIVI andfVl Eq. ([49]) is a quadratic in and its solutions are 

Q2 = 1 [Gn + G22 T \/a] , A = (Gii-G22)'+4Gi2G2i. (123) 

In what follows, the double valued odd order corrections l2n-i and S2n_i, n — 1,2,..., will be given for the upper 
sign in Eq. ([50]) . The sign ambiguity will refer to Eq. p23p . 
As a basis eigenvector we can take 

so(x) = {solix), so2{x)} - g{x) {1 , (Q^ _ Gn)/Gi2}, (124) 

where g{x) can be chosen in any convenient way. 

In hermitian vector cases, one should try to make ^1(2;) vanish by normalizing the above defined So(a;) and multi- 
plying it, if necessary, by the phase factor given by Eq. ([5T|) . i.e., by taking Sq = ei = e. 

In any of the theories under consideration, the basis vector in Ti.-^ is given by 

s^{x) = {-4,{x),s*,{x)}, (125) 

and 

si{x)^ci{x)s^{x), ci{x)^-2Q'D-^{s^,hm), (126) 
where D, given either by Eq. (|113p for hermitian cases or Eq. (|12ip in the non-hermitian theory, takes the form 

D = IsoiP G22 + |so2|' Gn - |soP - {s*^ so2 G12 + sqi s*02 G21) = ±\sq\^ VA. (127) 
Note that D is independent of the phase factor in Sg. 
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In Fulling's or Wronskian conserving hermitian theories we obtain (sq = e, s-*- = e-*-) 

s„(a;) = s,-^j(2;) + c,„(a;) e, (128) 

where Cm{x) = (e(x), Sm(a;)) is given by Eqs. p05p - (ll08p (with ei e). We reeah that Yo{x) = 1 and Yi{x) = 0. 
Note also that Q^D~^ is reaL 

In real hermitian cases, Sq, and G21 — G12 are all real, i.e., the asterisk can be dropped. 

In simplified theories (hermitian or non-hermitian) , Sm(x) is given by Eq. (|128p with Cm{x) = 0, where s-'-(x) is 
either normalized or non-normalized. 

In hermitian theories we obtain, see Eqs. ([70]) and ([74|) . 



Ym{x) = \so\-\so,h„^), (129) 
and in the non-hermitian theory, Eq. (|118p leads to 



1 



Q-^ci{x) G12IS0I2 



2soi 



In situations where we can eliminate the small parameter A by putting A = 1, |y„i(x)| should be small as compared 
to unity, see Eq. ([5^ . Furthermore, as in the lowest order s = Sg and Is-*-] = |so|, also the above defined multipliers 
c:^i{x) and Cm{x) should be small: 

|c;i;(a;)|, |c™(x)| «1. (131) 

If these requirements are not fulfilled, the phase integral approximation theory can only be used if we keep the small 
A parameter in the expansions of Y{x) and s{x). Examples of such behavior will be given in Sec. IVIIII 



VII. SINGULARITIES IN THE PIA AND CHOICE OF THE AUXILIARY FUNCTION a(x) 

For a given number of terms in the expansion pip , the applicability condition (|32p can always be fulfilled in any 
interval of x if the functions Ymix) are bounded there and we choose A sufficiently small. If Ym{x) are not only 
bounded but small as compared to unity, we can eliminate the small parmeter A by putting A = 1 . The analysis given 
in this section concerning the choice of auxiliary function a{x) is pertinent only to such cases. At the same time, the 
analysis concerning the singularities of the PIA is valid for any < A < 1. 

Putting A = 1 means that we treat both terms on the RHS in Eq. (|44)) or ((25)) on equal footing. In that case 
Eq. can be used as the definition of G(a:) or a{x) in terms of R(a;), i.e., for a given set of ODEs. 

Denote by Qo{x) an eigenvalue of the matrix R(a;) defining our set of equations ([3]). In the scalar case, Qo{x) = R{x). 
Eqs. (|44l) and (|48)) indicate that the matrices R(x) and G{x) have the same eigenvector s(x), and if A = 1, their 
eigenvalues Qq{x) and Q^{x) are related by 

Q\x) = Ql{x) - a{x). (132) 

In the scalar case, where only even order functions Y2n{x) are present, the question of accuracy of the PIA is 
addressed in Refs. 16 and 17. In Ref. 16, the behavior of Y2nix) in the vicinity of characteristic points is examined 
with particular reference to such points where Y2nix) — » 0. If that happens, the PIA tends to an exact solution 
of the wave equation ([1]). In Ref. 17, very accurate error estimates for the PIA are found, given in terms of the 
/i-integral introduced in Ref. 3, and two additional z/-integrals.^^. Generalization of these results to the vector case 
requires a separate treatment. However, as the scalar contributions to Ym{x) are always present in the vector case, 
the points where these contributions may be large or singular are also pertinent to the vector case. They correspond 
to singularities in cq (x) and can easily be determined by using simple models of the functions entering Eq. p32p . The 
power and exponential models discussed in Ref. 16 and 17 are both simple and useful. 

Locating the characteristic point for Q^{x) (zero or singularity) at a; = and denoting by Q'^{x) the model of 
Q'^{x) in some vicinity of a; = we can write 

Q^{x) = Qm{x) 1 + d{x) , d{x) ^0 as a; ^ 0. (133) 

If the small correcting function d{x) satisfies some additional (model dependent) requirements, one can relate a 
singular behavior of eg (x) at a; = to the behavior of (x) . Using then some transformation a; — s- i which maps the 
point X — into i = oo, we can also examine possible problems in the PIA at infinity. The simplest transformation 
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of this sort, x = l/cc, was shown in Ref. 16 to be the only transformation (up to translation and scaling) which 
can conserve ^^{x). However, it requires a simultaneous transformation of all functions entering Eq. (|132p . so as to 
conserve the products of times the function in question: 

x^Q^{x) = x^Q'^{5:), i.e., Q'^ii) ^ i-'^Q^{l/x}, (134) 

and similarly for Qq{x) and a(x), leading to eo(a^) = eo(S)- Using Eqs. ([18]), ((20l) . p6|) and (|133p we easily find 

SAQM]+a(x) d'{x)^lnQlix)-2d"{x) [d'{x)]' 



where "f{x) is a correction factor close to unity, 

j{x) = [I + d{x)]~\ (136) 

and Q^{x) is given by Eq. (|133p . This £9(2;) ^^id all three fractions in Eq. (|135p are invariant under the simultaneous 
transformation: 

{QmW: Qo{x),a{x)] ^ {Ql,{i), Qo{i)Mx)] - 3:-^{QU^/i), Qo{l/i),a{l/i)], d{x) ^ d(l/i). (137) 

In other words, this ea{x) and the fractions in question found in the limit a; — > are also valid in the limit x ^ 00 
and vice versa, if x and x are related by x a; = 1 . Eq. ()133p is also invariant under this transformation. 

As both X and d{x) are small, one can expect d'{x) and (i"(a::) to be simpler than d'{x) and d"{x). If that is the case, 
the last two fractions in Eq. (jl35p should be calculated from this equation rather than its equivalent for quantities 
with a tilde. An example of the opposite behavior is d{x) — exp(— a;) as i — > +cxd. 

The power model is defined by 

Q^(a:) = ca;™, c 0, i.e., Q^(5) = ci", m=-(m + 4), (138) 

«oo (^) = iq}^1+2 {"^("^ + 4) + 2 7(2^) ["^ X d'{x) ~ 2 d"{x)\ + 5 7' (a:) [x d'{x)\ ' } . (139) 

Here and in what follows, eoo(a;) = eo(a;)|a(a;)=o ^J^d the discussion will be performed in the x variable. It can be 
seen that the first term in braces, which represents the contribution to t^oix) coming from Q'^{x), is the leading 
contribution if m ^ 0, —4 and the derivatives of d{x) satisfy 

lim \x d'{x), x^ d"{x), . . .1 =0, equivalent to lim \x d'{x), x? d"{x), . . .1 =0. (140) 

X— 1-0 £^00 

In that case^^ the leading contributions to all functions Y2n{x) in the scalar case (not only to 12(2;) — eo{x)/2) come 
from Q'^{x), often also for a{x) 7^ 0. Note that Eqs. (|133p along with (|138p - (|140p are valid e.g., if Q^{x) is an analytic 
function which has a zero, a regular point or a pole at a: = . In that case, m is an integer, d{x) — X^fc^i '^kX^ 
and Eq. (|133p is a Taylor or Laurent expansion of Q^{x) about x = 0. In general, condition (|140p is fulfilled if the 
derivatives d'{x), d"{x), etc. are bounded as a; — > but it holds also for some singular behavior, e.g., for d{x) = a:lna;, 
^/x etc. Furthermore, the exponent m can be any real number. All these facts are also valid after the mapping x —^ x. 

The power model says nothing about eoo{x) if Q'^{x) or Q^{x) is an analytic function which has an essential 
singularity at x = or a; = c». This can be illustrated by the exponential models, given by 

Qliix) = cx~^ exp(77/a;), i.e., Quix) = c exp(77x), (141) 
eooix) = ~ 2 lix) [(4 x + v)x^ d'{x) + 2 x^ d"{x)\ + 5 ^Hx) [x^ d'(x)] (142) 

or 

Qliix) = c exp(7?/x), i.e., Qm{x) = cx~^ e:xjp{r]x), (143) 

^00 (x) = 5^^^\ , , Uir,-8x)-2j{x)[r,x^d'ix)+2x' d" (x)] + 5 7' (x) x [x^ d'{xf}, (144) 
loca:'^ exp(77/a;) ^. J 
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where ?] is a constant phase factor 
condition 



.]=0. 



(145) 



= 1. For these exponential models, condition (|140p can be replaced by a weaker 

\im\x'^ d'{x), x'^ d"{x), . . .] — 0, equivalent to lim \d'{x), d"{x),.. 

The RHS of Eq. (fT^ tends to zero as x if m + 2 < 0. That happens if 

lim x^ Q^{x) = lim x^ Q^{x) — oo, 

X — ^0 X — *C30 



i.e., if the limiting values of |(5^(a;)| or |(5^(i:)| are large as compared to \x\ 



(146) 

• This behavior is favorable for 

the PIA in the scalar case (I2T1 ^0). In the vector case, everything depends on the additional terms in Ym specific to 
this case, which will be illustrated by examples given in Sec. IVIIII For an analytic function, the condition m + 2 < 0, 
equivalent to to + 2 > 0, is fulfilled at higher order poles at x = 0, m = —3, —4, . . . (e.g., due to strongly singular 
potentials at a; = in the radial Schrodinger equation) or at simple zeros (to — —1), regular points (to = 0) and poles 
(to =1,2,...) at a; = 00. 

If the limits in Eq. (|146p are either finite or equal to zero, i.e., the limiting values of |(5^(a;)| or |(5^(i)| are not 



large as compared to \x\ 



or \x\ 



one can expect problems in the PIA both in the scalar and in the vector case. 



That happens for m + 2>0orm + 2<0. Thus, in the marginal case of to = 2 (equivalent to to = 2), eoo given by 
Eq. (|139p tends to — l/(4c), which can be large if |c| < 1/4. In the remaining cases of to + 2 > (i.e., 171 + 2 < 0), 
the RHS of Eq. (|139p tends to infinity, except for to = (i.e., ifi = —4). In that special case, in which Q{x) c ^ 
(and the corresponding Q^{x) has a zero of order four at i = 00), 



eooix) = ^{5j{x)[d'{x)]'~4d"{x)]}. 



(147) 



This is finite if the derivatives d'{x) and d"{x) are bounded (e.g., for analytic functions) but again can be large, 
especially if |c| is small. Singular behavior is also possible if the derivatives in question are singular. 

The relevance of condition (|146p to applicability of the PIA demonstrated here for the power model is in fact quite 
general. Its necessity, pertinent both to the scalar and the vector case, can be demonstrated as follows. If the limits 
in Eq. (|146p for Q'^{x) = Q'^{x) are either finite or equal to zero, so that condition p46p is violated, Q^{x) can be 
written 



Co + doix) , ^0(2;) 



as X 



where the constant cq may be equal to zero, leading to 

1 f 1 



eoo (a;) 



Co + do{x) { cq + do{x) 



[xd'{x)-x^ d"{x)] 



xd'{x) 



Co + (io(x) 



(148) 



(149) 



This tends to — l/(4co) if co ^ 0, i.e., can be large if |co| < 1/4, and is singular if co = 0. 

The sufficiency of condition (|146p for applicability of the PIA in the scalar case is conditional. Any function Qf^ix) 
for which this condition is fulfilled as x — > can be written 



Qn{x) — X b{x), lim b{x) 

x—>0 



00, 



leading to 



eoo 



ix) = -l[b{x)]-'\l + 



xb'{x) 


5 


\xb'{x)] 


2 

H 


x2 b"{x)\ 


b{x) 


4 


[ b{x) \ 




b{x) ] 



(150) 



(151) 



This eoo(x) will tend to zero if b~^ {x) x b' {x) / b{x) and b~^(x) x"^ b" {x) /b{x) tend to zero as a; ^ 0. That happens, 
e.g., for b{x) = x™, m < 0, and b{x) = exp(77/a;) ds, x ^ 0+, where 77 is a constant phase factor, \ri\ = 1, Rer^ > 
and TO is real. 

An important point is that if (5§(a;) has the form (|148p . the behavior of eoo(a^) given by Eq. (I149p . unfavorable for 
the unmodified PIA {a{x) = 0), can always be eliminated by an appropriate choice of the auxiliary function a{x). 
Thus, if we choose the leading term of a{x) to be proportional to x^^, 



a{x) = CaX ^[l+(ia(x)], Ca 7^ 0, da{x) ^ Q 



as X 







(152) 
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this term will also be the leading term of given by Eq. (|132p . This Q'^{x) will be given by Eqs. (I133P and (|138|) 

with 

TO = -2, C = CQ~Ca, d{x)^'^:^^ £at^a(^^ Ca ^ Cq. (153) 

Co - Ca 

Assuming that do (2;) and da{x) satisfy condition (jl40p . this condition will also be satisfied by d{x). In that case, using 
Eqs. (iMl), (fT39| and p52| we easily find 

eo{x) = |4c[l + d(a;)]} ^ [-1 + 4c^ - 7(2;) [a; d'(a;) + a;^ d"(a;)] +§7^(2;) [x d' (x)]^ + 4 Ca daix)] . (154) 
Choosing here and in Eq. (|152p Ca = 1/4, we obtain eo(x) — > 0, and Eq. (|132p gives 



^(2;)=Qg(a,)_ [i + d^a;)]^ i.e., g2(j) ^ g2(j) _ ^ ^ (^55^ 



For a finite zero or pole (at a; = 0), Eq. ()155p is applicable at zeros and the first and second order poles of Qq{x) 
(to > —2). For an infinite zero or pole {x — > 00), it can be used at higher order zeros of Qo{x) (m < —2). In either 
case, if TO, TO — —2, we must require that cq 7^ 1/4 in Eqs. (|133p and (|138p specialized to Q'^{x) — Q^{x). In the 
remaining situations, i.e., higher order poles at a; = (too < —2) and simple zeros, regular points or poles at a; = c» 
(toq > —2), the phase integral approximation with a{x) = ( "non- modified" approximation) in the scalar case tends 
to an exact solution of the wave equation. 

One should realize that in the case of a finite zero of Q^{x) (at a; = 0), the practical usefulness of Eq. p55p is 
limited by the fact that Q^{x) will have a simple zero in the close vicinity of a; = 0. This seldom happens in the 
remaining situations favorable for modification. 

For the Schrodinger equation in spherical coordinates, see Eq. ([7]) for = 1, the modification defined by Eq. (|155p 
with X = r and da{x) = gives justification to the well known "Langer modification";^ l{l + \) ^ {I + 1/2)^. 

In vector cases, apart from singularities in higher order corrections due to singularities in eo(x), as discussed above, 
additional ones occur at crossing points of the eigenvalues Q^{x). That is because elements of the A-"- matrix given 
by Eq. ([55]) contain the factor , where D is the determinant of the matrix — Q^l^ . This determinant vanishes 
at the crossing points, thereby introducing singularities in s,;^ given by Eq. (|68p . For N ^ 2, this is directly seen from 
Eq. (|127p . where A = corresponds to a double zero of the characteristic equation defining Q^{x). Denoting by Xcr 
the crossing point we obtain (A^ = 2) 

D{x)^g{x){x-x„T, 9{xcr)^0, (156) 

where p = 1. One should expect Eq. (|156p to be also valid if A^ > 2, but now with 1 < p < A^. 

With D{x) given by Eq. (|156p . s^^^ will contain the factor (a; — a;cr)~^ leading to a singularity at a; = Xcr- This factor 
will be small at sufficiently large distances from the singularity. At the same time the integral part of the coordinate 
(61,8™) given by Eqs. (|106p - (|108p wiU have a factor locally increasing with distance | ln|a; — a:cr| if p = 1. 

This, in contrast to scalar theory, may lead to situations where the higher order corrections are never small, unless 
we keep the small parameter A in the relevant equations. Note in this connection, that the integrals present in the 
remaining coordinates of Psm, (efc,Sm), fc > 1, will not contain the logarithmic contribution because the integrands 
in Eq. ([55)1 contain the derivative s:^'{x) rather than s;^(a;). 



VIII. EXAMPLES 



All examples in this section will be given for N = 2, where the algebraic part of the theory is simple, see Sec. I VII 
In the simplified theories (both hermitian and non-hermitian) which are algorithmic for A" = 2, all corrections in 
any order can be determined analytically. This may not be the case for FuUing's or Wronskian conserving theory 
due to the integrations present in Eqs. (|106p ~ p08p . Nevertheless, one can write a symbolic program in Mathematica 
pertaining to all theories^" which can be useful in many cases. All results presented in what follows were produced 
by using this program. In all examples, we choose a{x) = and the eigenvalues are real. 
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A. Simple real hermitian case with Q2 > 



The eigenvalue given by Eq. p23p takes a simple form if A is a perfect square. An example of this type (real 
hermitian case) is^^ 



„ , > , a; cos X + sm x (x — 1) cosx sma; 
' (a; — 1) cosx sma; a; sm a; + cos a; 



A = (a; -1)2, 



(157) 



where Gn + 6*22 = a; + 1, and so either = 1 or = x. 

By slightly modifying this G matrix we can produce simple examples of other vector cases (real hermitian with 
< 0, complex hermitian and non-hermitian) which will be discussed in the following subsections. Here and in the 
following simple examples related to Eq. (|157p we assume a; > 0. Note that the eigenvalues are much simpler than 
the elements of the G matrix, which is rather unusual (no complication due to the square root in Eq. (|123p ). Another 
peculiarity of this example is that all integrals in Eqs. (|106p - (|108p . as well as those in Eq. can be calculated 

analytically. Using the program in Mathematica based on Eqs. (|123p - (|130p .2° we can determine analytically all 
quantities pertaining to any (reasonable) order of both FuUing's and Wronskian conserving theory. We recall that in 
all hermitian theories, the basis can be chosen so that Yi{x) = 0. 

As the eigenvalue = 1 is a; independent, i.e., eo{x) = 0, all scalar contributions to l2n(a;) are zero. That 
means that only terms specific to the vector case in the higher order corrections Ym{x) are present. Both types of 
contributions are present for — x. In that case, we are dealing with a simple zero of Q'^{x) at a; = and the first 
order pole at x = 00. The function eo(x)(= 5/(16x'^)) has a pole at x = (unfavorable behavior at a zero of Q'^{x)) 
and tends to zero as x — > 00 (favorable behavior at a pole of Q'^{x)). For both eigenvalues, a singularity should occur 
at their crossing point, x = 1. 

In FuUing's hermitian theory, Y2n{x) are real as expected and Y2n-i(x) = as required for current conservation. 
(We recall that the FuUing's hermitian theory requires reality of all corrections Yjn{x), whereas in a real hermitian 
case, Y2n-iix) are pure imaginary.) 

If = 1 we obtain 



e = {sinx, — cosx}, = {cosx, sinx}, £0(2;) = 0, 



(158) 



Yi{x) = 0, 



X - 1 



ciix) 
c^(x) 



2i 
~ x-1' 
4 

{x-iy 



8 ln|x - 1| 
X- 1 



Ci(x) = — 4i ln|x — 1|, 



(x-iy 



Inix- 11 



(159) 
(160) 



Far away from the singularity (x 3> 1), only the multipliers c:^{x) are small and tend to zero as x — > cx) (the logarithmic 
terms are always divided by positive powers of (x — 1)). In the same limit, the corrections |y2Ti(a;)| are smaller than 
unity and decrease with increasing n but tend to non- vanishing constants {— 5, g , yfg, • ■ •) and the multipliers 
c,„(x) are logarithmically divergent. This means that the higher order corrections in s(x) for x 3> 1 can only be small 
(in any finite interval of x) if we keep A in the expansions pip and (|47p and take it sufficiently small. This seems to 
be quite a general feature of FuUing's theory if Q^{x) > 0, see e.g., the results below. 
If = X we obtain: 



e = {cosx, sinx}, e ={— sinx, cosx}, eo(x) 



16x3' 



(161) 



Y,{x) 

C2{x) 



0, c^{x) 



2U 



x-1' 
2 5 1 

~ 2x 



ci(x) = 4i 2^/x ~ h(x) , h{x) — In 
32x'' + 64x3 



^/x + l 



x-1 32 x3 
-16x2 + 32 X- 17 



2x- 



{x-lf 



2x (x 



8/i(x) h{x) - 4v^ 



29x2 + 6x-l S^x 
^ X - 1 



1 



ft,(x). 



(162) 



(163) 



The first and third term in Y2{x) are specific to the vector theory, whereas the second term {— eo(x)/2) is the scalar 
contribution. All tend to zero as x 00, and this is true for contributions to any Y2n{x). As for the velocity 
multipliers, only c^{x) — s- as x — s- 00. In the same limit, c^(x) — s- c 7^ 0, |c^(x)| — > 00 for n > 2 and |cm(x)| — s- 00 
for any n. 
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In the Wronskian conserving theory, l2n(a;) are again real, Yi{x) = and there are no integral contributions to 
(e, s^), as expected. However, the requirement 1271-1(2;) = is violated if n > 1. Therefore, the Wronskian for u*(.t) 
will be conserved through second order only. The relevant results are as follows. 

If = 1 we obtain 

Yi{x) = 0, ci{x) = ci(x)=0, (164) 



4 



Y2{x) = -TTTI^' c^(^) = 7—^' C2(x) - 7—^, (165) 



Yiix) = — 7T3-,---, (166) 

and the results for = x are: 

Yi{x) = 0, dix)^^, ci(a;) = 0, (167) 
X ~ 1 

... N 2 5 1 I , , 3x^ + 6x-l , , 2x 

^^^"^ = ^+32^-2^' ^^("^= 2.(.x-l)3 ' ^^(")=(^^' (1^^) 

In the simplified hcrmitian theory, the corrections through 2nd order are obtained if we put ci{x) = in the above 
results for the Wronskian conserving theory. Next order corrections have a similar structure in both theories but are 
not identical. The odd order corrections, Ys, ¥5, etc. are non-zero, e.g., 

r3(x) = -2^^^±lJ, (170) 
(a; — 1)-^ 

if = 1, etc. They are all pure imaginary, as expected, and so after integration (i J Y2n-i{x) Q{x) dx) will contribute 
to the amplitude rather than the phase of u^(x) given by Eq. (pT|) . The behavior of higher order corrections is now 
fully analogous to that in the scalar case (tend to zero as x ^ 00) except for the (somewhat peculiar) behavior of 
Y2n{x) if = 1, the same as that in FuUing's or Wronskian conserving theory. The basis vectors e and e-'- arc the 
same in all theories. 



B. Simple real hermitian case with < 
By changing sign of the diagonal elements in the G matrix (|157p . i.e., for 




(171) 



the eigenvalues change sign, i.e., either = — 1 or = —x. The corrections for the resulting real hermitian case with 
< are very closely related to those for the original real hermitian G matrix with > 0. They are all real and, 
except for differences in signs and the absence of i for < 0, the corrections in the Wronskian or current conserving 
theory for < are the same as those in the current or Wronskian conserving theory for > 0, respectively. And 
in the simplified hermitian theory, the corrections in both cases are the same in this sense. However, after dividing 
the odd order corrections for > by i, one has to replace 

{ei, 62} {61,-62} and {6]^, e^} {-6]^, 6^} 

and change signs in y,„ and c„i for m mod 4 = 1 or 2 and in c^, for m mod 4 = 3 or 0. In particular this means that 
the required and actually observed vanishing of y2n-i(a;), n = 1, 2, . . ., in FuUing's theory for Q^{x) > implies the 
analogous vanishing of Y2n-i{x) in the Wronskian conserving theory for Q^{x) < 0. In other words, there should be 
no problem in using the Wronskian conserving theory for Q^{x) < in contrast to Q^{x) > 0, where it can be used 
through the second order only. 

If = — 1, the results of the Wronskian conserving theory are, in view of Eqs. (|158P " (|160p . 



e = {sinx, cos x}, e ={— cos x, since}, 



(172) 
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Yi{x) = 0, ciix)^ -, ci(x) =4 1n|2;-l|, (173) 

a; — 1 



= c.(.) = ^+8 1n^|.-l|,.... (174) 



Similarly, the results for ~ ~x will be closely related to those given by Eqs. p6ip ~ (|163p . etc. 

C. Simple complex hermitian cases 

For any hermitian matrix G(x), on multiplying 6*12(2;) by a constant phase factor e*''' and 021(2;) by e~^^ , where 
tf is real, another hermitian G matrix is obtained. With this transformation, the eigenvalues Q^, the corrections 
Yra{x) in all theories and the multipliers Cm = (e, s^) given by Eqs. (|106p - (|108p are left unchanged. In the remaining 
results, simple replacements only are needed: 

{61,62} ^ {e"^ 61,62}, {6]^,6^} ^ {e]^,e"*'^6^}, c:^^ e'"^ c^. 

They leave |e| and le-*"! unchanged. 

If one starts with the real hermitian 0(2;) matrix (|157p or p7ip . one obtains simple examples of complex hermitian 
cases (e.g., for e^'^ = i or (1 + i)/\/2, etc.). Unfortunately, they are not "complex enough" so as to avoid making 
the RHS of Eq. (|94|) for < equal to zero identically in the simplified hermitian theory (which happens if 3(2;) is 
real). Non-trivial examples require if to he x dependent. For the simplest choice, (p = x, the integrand in Eq. (j8ip 
is non-zero but the integral is elementary. This means that the higher order corrections in the simplified hermitian 
theories (with 11(2;) = 0) can easily be determined. Unfortunately, the integrals in either the current or the Wronskian 
conserving theory, Eqs. pOSp - pOSp . are non-elementary. 



D. Simple non- hermitian case with Q2 > 

On multiplying 6*12(2;) of the G matrix (jl57p by any real or pure imaginary number or function and at the same 
time dividing 621(2;) by the same quantity, the eigenvalues Q^{x) and A given by Eq. p23p will be left unchanged. 
An example of the resulting non-hermitian matrix is 

„, ^ ( X cos^ 2; + sin^ X 2 i (2; — 1) cosx sinx\ 

G(2;) = .1/ ix • -2 9 ■ (175) 

I ^ (2; — 1) cosx sma; 2; sm 2; -I- cos 2; ) 

If = 1, it is convenient to choose g(x) — 2 sinx in Eq. (|124p . as it eliminates the singularity in the basis vectors 
and leads to Yxix) = (rather unusual in a non-hermitian case). With this choice we obtain: 

So = {2 sinx, i cosx}, s-'" = {i cosx, 2 sinx}, £0(2;) = 0. (176) 

The results of our non-hermitian theory are: 

ri(x) = 0, ci{x)^— — , d(x) = 5-3cos2x, (177) 

(x — 1) d(x) 

Y^{x) = — \^ , [-59x^-f 26x + 33 + 60(x- l)^cos2x 

4 (x — l)^a^(x) 

-9(x2 + 2x - 3) cos 4x + 12(10 sin 2x - 3 sin 4x)] , 

If = X, a convenient choice in Eq. (|124p is ^(x) = 2 cosx, again eliminating the singularity in the basis vectors 
and leading to yi(x) = 0. With this choice we obtain: 

5 

So = {2 cosx, —i sinx}, s^ = {— i sinx, 2 cosx}, £0(2;) = , (179) 

lox'^ 
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and the results of the non-hermitian theory are: 

Yi{x) = 0, CjL(x) = , ^i^u ^ ^ d{x) = 5 + 3 cos2x, (180) 
(a; — 1) d(x) 

Y2(x) = — — ---[5282;'' + 416x3 - 6492:2 - 590j: + 295 -(960x^- 1920 

64 x"^ [x — 1)^ a^(xj 

+660 + 600 X - 300) cos 2x + (432 x^ - 288 x^ - 99 x^ 
-90 X + 45) cos 4x + x^ (x + 1 ) (960 sin 2x + 288 sin 4x)] , 

c^(x) = — ^ ^3 [- 15 x^ - 30 X + 5 - 3 (3 x^ + 6 X - 1) cos 2x + 24 x^ (x - 1) sin 2x] , . . . (181) 

Next order corrections can be generated by using an appropriate program in Mathcmatica.^*' AU are oscillating and if 
= X, the amplitude of the oscillations tends to zero as x ^ oo. Therefore in that case, all higher order corrections 
are small for x ^ 1. 

Note that in the previous examples pertaining to hermitian cases, only the basis vectors contained the sin and cos 
functions and were therefore oscillating like the the corresponding G matrices. No such functions were present in 
Y2n{x) or the multipliers c;^(x) and c„(x). In view of Eqs. p27p and (|129p . an essential point was that the basis 
vectors e and e-"- were normalized, in contrast to those given by Eqs. p76p and (|179p . 

E. More realistic real hermitian case with < 

When numerically examining small oscillations of a single quantum vortex in Bose-Einstein condensate jiSiii^ one 
arrives at the eigenvalue problem for a set of two one-dimensional Schrodinger like differential equations: 



d u 1 du 

dr^ r dr 

d^v 1 dv 

dr'^ r dr 



2(j)l{r) + ^^l + 2uj + k^ u- 4>l{r)v = 0, 
[202 (r) - I - 2uj + k'^]v ^ <t)l{r)u = 0, 



(182) 



in which r is the cylindrical radius (0 < r < oo) k is the wavenumber and m is the frequency of the oscillations 
{k, uj > 0), and (j)o{r) is the radial profile of the vortex. This profile is a monotonic function described by a nonlinear 
2nd order differential equation and subject to the boundary conditions 0o(O) ~ 0, and 4>{r) — > 1 as r — *■ (X). To find 
initial conditions for numerical integration of the differential equations p82p at some large but finite value r = ras, 
these equations were first transformed to the reduced form ^ in which (wi(x) = x^/'^u{x), U2{x) — x^^'^v{x), x — r) 

^ ' \ h2(x) ho(x) + hi{x) J ' ^ ' 

ho{x) = -l-k^ + doix), hi{x) = 2{LU + x-^), h2ix) = -I + di{x), (184) 

^ , , 1 4 38 748 , , , 1 2 19 374 . „ , 

doix)^^ + ^ + ^ + -^, di{x)^^ + ^ + -^ + ^, (185) 

where Eqs. (|185p follow from the asymptotic expansion of (j)o{r) as r — > oo. Choosing a(x) = and A = 1 (allowed 
in the simplified hermitian theory), the eigenvalues of the matrix G(x) = R(x) given by Eq. (|183p can be written 
iQ' < 0): 

\Q\ = ^-hoix)±r{x), r(x) = ^hl{x) + hl{x), (186) 
and the corresponding eigenvectors So(x) are 

/ \ / \ f-, hi(x)Tr(x)) /-,„„N 
So(x) =ff(x) 1, ' y (187) 



h2{x) 

where the factor g(x) can be chosen in any convenient way. In Refs. 10 and 19, .g(x) = 1, leading to 

''^(^)=^^^(^)2^^(§^5^M^' ^^i^)^h,{x)h'2{x)-h2{x)h[{x). (188) 
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Second order corrections are much more complicated. 

It can easily be shown that cj-{x) is independent of the factor g{x) in Eq. (|187p and thus is given by Eq. (|188p for 
any g{x). In contrast, Yi{x) depends on the choice of g{x). In particular we obtain Yi{x) = if So{x) is normalized. 

In Refs. 10 and 19, an approximate analytical solution of Eq. ^ tending to zero as a; — *■ cxd was looked for in the 
form of a linear combination of two zero order phase integral approximations u~(a:), see Eq. (j6ip {Q'^{x) < and 
A = 1). These approximations were referred to as Ugo ("greater exponent", for the upper sign in Eq. ()186p with 
hQ{x) < 0) and Ugc ("smaller exponent"): 

u(a;) = Csc Uso(x) + Cgc Ugc(a;). (189) 

Eq. (jl89p was used for x > Xas, with Xas determined experimentally as the minimal value above which the eigenvalue 
uj{k) was practically independent of Xas- We were interested in uj{k) in the k limit, which in general required 
a^as > 2.2/ k. Here, by going to higher orders of the simplified hermitian theory, we can determine quantities related to 
the error of the lowest order approximation. The relatively simple eigenvectors (|187p with g{x) = 1 are not normalized 
which makes Yi{x) non vanishing. Normalization of So(a;) introduces a complication in the lowest order but simplifies 
formulas for higher order corrections (due to Yi{x) = 0). Using our present results, we can compare numerical values 
of the first and 2nd order corrections in Y{x) and s(x) at x = Xas for the normalized and non-normalized eigenvectors. 
These corrections tend to zero as x — > c» and can be expected to fall below their values ai x — x^s. That is because 
the R matrix (|183p tends to a constant matrix as x ^ oo. In particular this implies 



lim \Q{x)\ = \/l + A:2 ± ^1 + 4cj2 ^ J y2 ^ Um so(x) = {1, -2w ± Vl + 4^2} ~ {1, ±1}, (190) 

where the RHSs give the leading term in the fc — > limit in whic h^^i^^ cu ~ i/c^ ln(l/A:) also tends to zero. The limiting 
value of |(3(x)| for the lower sign is small which is unfavorable for the PIA (as it makes eo large, see Table |T|. This 
was the actual source of problems in the numerical solving of the eigenvalue problem in question. The minimal value 



TABLE I: First and 2nd order corrections for non-normalized and normalized eigenvectors at s = 55, k = 0.04 and uj = 0.002604. 
The first two lines represent Ugo and the last two lines Use. 



\Q\ 


eo/2 


Yi 


Y2 


ci 




1.41464" 


-2.54639-10"* 


-8.4752-10"'' 


1.3783-10"'^ 


-1.70658-10"^ 


-9.88846-10"^ 


1.41464'' 


-2.54639-10"* 





-2.55731-10"* 


1.70658-10"'' 


-9.88846-10"^^ 


0.0427842" 


1.59832-10"^ 


2.83539-10"'' 


1.58104-10"^ 


5.16137-10"'^ 


-3.15819-10"'^ 


0.0427842' 


1.59832-10"^ 





1.59832-10"^ 


5.16137-10"^^ 


-3.15819-10"'^ 



"Non-normalized eigenvector. 
'Normalized eigenvector. 



of k which could be treated in our calculation was k — 0.04. The corresponding numerically found quantities were 
u! = 0.002604 and Csc/Cgc — —4.367. To get uj practically independent of Xas required Xas = 55. For these values, the 
contribution to u(xas) coming from the first term in Eq. p89p is nearly 26 times larger than from the second term 

(|usc(xas)|/|ugo(xas)| — \J V2 / k = 5.946). One can easily see that this fact also holds for the normalized eigenvector. 
This means that the relative error in u(xas) is defined by the relative error in Usc(a;as)- This in turn is given by the 
corrections Yi(a;as), Y2{xas) and the multipliers c;j^(xas) and c^(a;as), see the last two lines in Table ID The dominating 
quantity is l2(2;as)- It depends very little on normalization. Therefore one should not expect an improvement in the 
zero order approximation due to normalization of So(x). 

Note that in most cases, the second order corrections l2(a;as) are very close to their values in the scalar case, 
eo(a:as)/2. 

IX. CONCLUSIONS 

We present four generalizations of the Phase Integral Approximation to sets of ODEs of the Schrodinger type (three 
for hermitian and one for the non- hermitian sets). 

The first is an extension of FuUing's hermitian theory^ which conserves the generalized current in each expansion 
order. In Ref. 8, this theory was developed for positive definite matrices. In that case and for real hermitian matrices, 
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this theory conserves both the current associated with the complex PIA and the Wronskian built from its real and 
imaginary part. Using it for negative definite matrices makes the current vanish (and thus be conserved). This can 
be of some interest only in such complex hermitian cases with < in which in the simplified hermitian theory, 
the current given by Eq. (|94[) is not identically zero. This current is identically zero in real hermitian cases where the 
s(x) vector is real, but also in simple complex cases, see Sec. IVIirCl for examples. 

The second generalization is the hermitian theory that conserves the Wronskian built from two PIAs u^{x). Our 
examples show that applicability conditions for this theory are fulfilled through second order only if > 0. No 
such restrictions were found for negative definite matrices. In that case and for real hermitian matrices, this theory 
conserves both the Wronskian built from (real) u^(x) and the current associated with an approximate complex 
solution u+(x) + iu~{x). 

The third theory ( "simplified hermitian" ) conserves the current and the Wronskian only in lowest order but contains 
no integrations characteristic of the first two. In the non-degenerate case, it contains no integrations in higher order 
corrections and is thus fully algorithmic. Furthermore, these higher order corrections were never found to be large far 
away from singularities. In such a situation one can eliminate the small parameter A by putting A = 1 and the idea of 
modification described in Sec. IVIll is applicable. This theory is the simplest and in applications (like that described 
in Sec. IVIII E|l may turn out to be the best. For N = 2, using the program in Mathematica described in Ref. 20, one 
can determine corrections of any (reasonable) order. No matter how complicated they are, their numerical values can 
be determined. Writing a similar program for iV = 3 is not difficult. For more than three equations, linear algebra 
becomes more complicated and finding the eigenvalue analytically may be impossible. 

The non-hermitian theory for non-degenerate vector cases developed in Sec. |V] is based on the same assumption 
as that in the simplified hermitian theory ((ei,Sm) = 0). The example given in Sec lVIII Dl seems to suggest that 
in typical situations {Q^ = x), higher order corrections are small at large distances from singularities {x —> oo). 
This theory is fully algorithmic and the idea of modification should be applicable. Obviously, it can be used also 
for hermitian non-degenerate matrices. In that case it gives the same results as the simplified hermitian theory but 
without the necessity to introduce the basis in the orthogonal complement of the eigenvector. 

This paper only deals with the adiabatic part of the PIA for the vector case. The connection problem^i^iiii^ requires 
a separate treatment. In the references just mentioned, this problem for the scalar case was solved by tracing the 
unknown function u{x) in the complex x plane. Note in this connection, that our non-hcrmitian theory developed in 
Sec. |V] is the only vector theory of the PIA that allows for complex x. Only within this theory can one try to solve 
the connection problem for the vector case by tracing the unknown vector u(a;) in the complex x plane. 
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